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This  paper  presents  an  algebraic  method  for  parameter  identification  of  Thevenin's  equivalent  circuit 
models  for  batteries  under  non-zero  initial  condition.  In  traditional  methods,  it  was  assumed  that  all 
capacitor  voltages  have  zero  initial  conditions  at  the  beginning  of  each  charging/discharging  test.  This 
would  require  a  long  rest  time  between  two  tests,  leading  to  very  lengthy  tests  for  a  charging/discharging 
cycle.  In  this  paper,  we  propose  an  algebraic  method  which  can  extract  the  circuit  parameters  together 
with  initial  conditions.  This  would  theoretically  reduce  the  rest  time  to  0  and  substantially  accelerate  the 
testing  cycles. 

©  2014  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Several  types  of  battery  models  are  used  to  evaluate  the  per¬ 
formances  of  batteries  and  to  study  their  interactions  with  other 
electrical  devices.  As  compared  with  electrochemical  and  mathe¬ 
matical  models,  the  electric  circuit  based  models  are  much  simpler 
for  computation  and  analysis.  They  have  a  high  potential  in  terms  of 
accuracy  and  parameterization  effort  [1],  They  are  the  most  intui¬ 
tive  for  use  in  circuit  simulations  [2]  and  can  be  easily  integrated 
with  other  electric  devices  and  control  algorithms  at  a  system  level 
[6], 

A  widely  used  circuit  based  models  for  batteries  are  Thevenin’s 
equivalent  models  as  depicted  in  Fig.  1. 
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It  was  initially  proposed  in  Ref.  [24]  in  1992  for  led-acid  batteries 
and  then  widely  used  in  Refs.  [8—34]  to  model  various  types  of 
batteries  such  as  lead-acid,  lithium-ion  (Li-ion),  Li-polymer,  nickel 
metal  hydride  (NiMH),  and  fuel  cells. 

Thevenin's  equivalent  models  have  been  used  for  analysis, 
design,  and  simulation  of  battery  powered  electronic  systems 
[8,10,11,16—18,20,33].  They  are  also  important  for  characterization 
of  battery  performance,  life-time  estimation,  power  management, 
and  efficient  use  of  batteries  [21—25],  In  Ref.  1],  three  Thevenin's 
equivalent  circuit  models  were  applied  to  simulate  a  real-life 
driving  cycle  of  an  electric  vehicle.  In  Ref.  [4],  the  models  were 
used  for  online  estimation  of  state-of-charge  and  open  circuit 
voltage  of  lithium-ion  batteries  in  electric  vehicles. 

The  circuit  model  in  Fig.  1  has  several  variations.  For  example, 
the  ideal  voltage  source  can  be  replaced  with  a  capacitor  for 
describing  the  charging  dynamics  or  longer  time  discharging 
behavior.  The  ideal  voltage  source  can  also  be  replaced  with 
another  pair  of  capacitor  and  resistor  in  parallel  to  account  for  the 
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Fig.  1.  Thevenin's  equivalent  circuit  model  for  batteries. 


Fig.  2.  A  battery  connected  to  an  electronic  load. 


self-discharging  behavior.  These  variations  have  also  been  widely 
considered  in  the  literature  (e.g.,  see  Refs.  [9,24]).  More  advanced 
circuit  models  have  been  developed  based  on  Fig.  1.  For  example,  in 
Ref.  [35],  the  independent  voltage  source  is  replaced  by  a  controlled 
voltage  source  which  is  dependent  on  the  voltage  of  a  capacitor  in 
parallel  with  a  resistor  and  a  dependent  current  source. 

It  is  known  that  the  parameters  in  the  circuit  model  depend  on 
many  factors,  such  as  state  of  charge  (SOC),  the  load  current,  the 
temperature  and  even  the  history  of  charge  and  discharge 
[8,20,24,26].  They  can  be  regarded  as  constants  under  a  certain 
working  condition  and  over  a  relatively  short  period  of  time.  More 
complex  nonlinear  models  have  been  proposed  based  on  the  rela¬ 
tionship  between  the  circuit  parameters  and  the  operating  condi¬ 
tions  [11,14,15,20],  In  Ref.  [5  ,  the  resistances  Rg  and  R2  were 
dependent  on  the  open  circuit  voltage  and  the  capacitor  voltage. 

A  core  task  in  battery  modeling  is  to  identify  the  parameters  in 
the  Thevenin's  equivalent  circuit  model  via  some  experimental 
responses  under  a  given  operating  condition.  With  families  of  pa¬ 
rameters  identified  for  various  working  conditions,  a  nonlinear 
model  describing  the  dependence  of  these  parameters  on  SOC, 
current,  and  temperature,  can  be  obtained  with  numerical  methods 
[14],  In  Refs.  [2,5],  multi-objective  genetic  algorithms  were  devel¬ 
oped  for  extracting  the  parameters  of  the  linear/nonlinear  circuit 
models. 

The  circuit  parameters  are  traditionally  identified  by  using  least- 
square  curve  fitting  of  the  experiment  data.  In  Refs.  [12,13],  we 
developed  simple  analytical  methods  for  identifying  the  parame¬ 
ters  for  several  equivalent  circuit  models  including  the  one  in  Fig.  1. 

A  complete  nonlinear  or  parameter  dependent  circuit  model  for 
a  battery  requires  many  rounds  of  charging  or  discharging  cycles. 
Each  round  consists  of  many  tests  for  the  characterization  of  the 
parameters'  dependence  on  SOC.  Between  two  subsequent  tests,  a 
sufficient  rest  time  is  required  due  to  the  assumption  of  zero  initial 
condition  at  the  beginning  of  each  tests.  The  rest  time  could  be 
30  min  or  longer,  depending  on  the  type  of  battery  [36], 

In  this  paper,  we  propose  an  algebraic  method  for  extracting  the 
circuit  parameters  together  with  the  initial  conditions.  The 
assumption  of  non-zero  initial  condition  will  eliminate  the  need  for 
rest  time  or  reduce  the  rest  time,  thereby  accelerate  the  modeling 
process. 

2.  Algebraic  method  for  extracting  the  parameters  and  the 
initial  conditions 

We  will  first  show  that  it  is  impossible  to  extract  the  parameters 
and  the  initial  conditions  at  the  same  time  when  a  constant  load 
current  is  applied.  If  the  load  current  is  a  step  function,  i.e.,  it  takes 
one  value  over  a  period  of  time  and  then  switch  to  another  value, 
then  the  parameters  and  the  initial  conditions  can  be  identified 
simultaneously. 


2.3.  Parameter  identification  via  a  constant  load  current 

For  simplicity  and  without  loss  of  generality,  we  consider  the 
battery  model  with  2  pairs  of  parallel  resistors  and  capacitors, 
(Ri,Ci),(R2,C2),  as  depicted  in  Fig.  1.  This  will  be  called  the  2nd-order 
model.  The  method  can  be  easily  extended  to  3rd  or  higher  order  by 
using  the  algebraic  method  in  Ref.  [12], 

A  common  setup  to  identify  the  parameters  is  to  connect  the 
battery  to  an  electronic  load  which  absorbs  a  prescribed  current  i 
from  the  battery,  see  Fig.  2. 

The  current  can  be  a  constant,  a  step  function,  a  square  wave,  or 
other  types  of  programmable  functions.  The  terminal  voltage  v(t)  is 
recorded  and  used  for  parameter  identification. 

In  this  section,  we  consider  a  constant  load  current  i  =  I  and 
assume  that  the  load  is  connected  at  t  =  0  (implying  i  =  0  for  t  <  0). 
A  traditional  assumption  made  on  the  circuit  is  that  the  initial 
conditions  vi(0),v2(0),...,  are  all  0.  In  this  paper,  we  don't  make  this 
assumption  and  allow  the  initial  conditions  to  be  nonzero  and 
unknown.  Denote 

vt  (0)  =  fio,  f2(0)  =  i’2o. 

The  instantaneous  terminal  voltage  before  applying  the  current 
load  is 

v(0-)  =  E  -  vw  -  v20  (1) 

Under  a  constant  current  i  =  I,  the  terminal  voltage  response  for 
t  >  0  is 

v(t)  =  E  -  R0I  -  R^I  -  R2I  -  (v10  -  R]l)e  ^  -  (v2o  -  R2l)e^v^ 

(2) 

Comparing  (1)  and  (2),  we  have  Rg  =  (v((T)-v(0+))//. 

The  other  parameters,  £,  Ri,R2,Ci,C2,  and  initial  condition  Vi0,v2o 
need  to  be  computed  from  the  voltage  response  v(t).  Recall  that 
under  the  assumption  of  zero  initial  condition  v10  =  0,v2o  =  0,  all 
the  parameters  can  be  identified  from  the  terminal  voltage 
response  using  various  methods.  However,  when  vio,v2o  are 
nonzero  and  also  unknown,  the  parameters  cannot  be  identified 
based  on  the  response  in  (2),  in  spite  of  obtaining  infinitely  many 
equations  for  only  7  unknowns.  This  is  because  the  parameters 
cannot  be  separated  from  the  initial  condition  under  a  constant 
discharging  current.  This  fact  is  formally  stated  in  the  following 
claim. 

Claim  1.  Let  the  terminal  voltage  v(f)  be  given  for  t  >  0.  The  var¬ 
iables  E,  R1,R2,Ci,C2,Vio,v2o  cannot  be  uniquely  determined  from  v(t). 
Denote  si  =  vio-Rih  s2  =  v2q-R23  and  t  1  =  R1C1,  r2  =  R2C2. 
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1.  S\,S2J\J2  can  be  uniquely  determined  from  5  points  of  v(t). 

2.  For  the  5  unknowns,  E,Ri,R2,Vio,v2o,  there  are  only  three  linearly 
independent  equations. 

A  constructive  proof  for  the  claim  is  given  below.  We  first  show 
item  1.  Denote  si  =  vio-Ril,  S2  =  v2o-R2E  Choose  T  >  0  and  let 
dj  =  e~T/R,Cl , d2  =  e~T/R2c2.  Then  for  any  integer  k> 0, 

v(kT)  =  E  —  R0I  —  R^I  -  R2I  —  Sjd|  -  s2d\, 

The  constant  term  E-R0I-R1I-R2I  can  be  eliminated  by  sub¬ 
tracting  v((k  +  1  )T)  from  v(kT).  To  be  specific, 

v(o+)-v(r)=s1(d1  -i)  +  s2(d2-i) 
v(T)  -  v(2T)  =  s1(d1  -  l)dt  +  s2(d2  -  1  )d2 


v(2 T)  -  v(3T)  =  Sl  (d,  -  l)d?  +  s2(d2  -  1  )d\ 


2.2.  Parameter  identification  via  a  step  load  current 


A  simple  approach  to  introduce  more  linearly  independent 
equations  is  to  change  the  value  of  the  current  to  a  different  con¬ 
stant.  To  be  specific,  we  choose  a  switching  time  T0  and  let 


fit)  = 


1 ,  0  <t<T0 

h,  t>To 


We  can  use  the  measurement  from  0  to  To  to  find  Si,s2,ti,t2  and 
use  the  measurement  after  To  to  form  two  more  linearly  indepen¬ 
dent  equations. 

Claim  2.  All  the  circuit  parameters  and  the  initial  conditions  can 
be  identified  from  measurement  taken  by  applying  a  step  load 
current. 

A  constructive  proof  for  this  claim  is  given  below. 

Choose  T  <  To/4.  We  can  use  the  method  in  the  previous  section 
to  compute  si  =  Vw-Ril,  s2  =  v2o-R2I  and  ri  =  R1C1,  t2  =  R2C2  by 
using  v(f)  for  t  <  To. 

Denote  aq  =  e~r°/n ,  a2  —  e~T°/T2.  Then  at  To,  the  two  capacitor 
voltages  are 


v(3T)  -  v(4 T)  =  s,(d,  -  1  )df  +  s2(d2  -  1  )d\ 

Define  x\  =  Si(di-l),  x2  =  s2(d2-l),  the  above  equations  can  be 
rewritten  as 

*i  +x2  =  i/(0+)  -  v(T) 
x1d1+x2d2  =  v(T)-i/(2T) 
xxd\  +x2dl  =  v(2T)-v(3T) 


vi  (T0)  =  s^a  1  +  R^I, 

v2(T0)  =  s2a2  +  R2I. 

Choose  Ti  >  0.  Let  p1=e~Tl/T>,  p2  =  e^T'lT1.  Then  for 

k=  0,1,2,..., 

vi (T0  +  kh)  =  (tqOTo)  -  R,h)p\  +  R,I 

=  (s1«1+R1(/-/1))p'[+R1/ 

v2(Tq  +  kh)  =  (v2(T0)  -  R2fi)pk2  +  R2I 

=  (s2a2  +R2(I  -  h))p2  +  R2I 

Denote 


x^d\  +x2d\  =  v(3 T)  -  i<(4T) 


Y\  =si«i  +R1(I-h) 
y2  =  s2a2  +  R2(I  -  I\). 


From  the  above  equations,  the  4  variables  x^x2,d\,d2  can  be 
easily  solved  by  using  the  algebraic  method  in  [12]  (details  will  be 
provided  in  Section  2.3).  Then  Si,S2,ti,t2  can  be  computed  from 
these  variables: 


Si  = 


*i 

di  — r 


S2 


*2 


d2-r 


We  have 

v(T0  +  kTi)  =  E  —  Rofi  -  Rih  -  R2h  -yiPi  ~y2p2 
By  subtracting  two  subsequent  points,  we  obtain 

v(T+)  -  v(T0  +  Ti)  =yi(pi  -  1)  +y2(P2  -  1) 


(4) 


ri=*lCl=-h5?  T2=Rl°2  =  ~hk' 

With  Si,s2,T1,r2  uniquely  determined,  we  have 
E  -  Ri  I  -  R2I  =  v(t)  +  RqI  +  s1e~t/Tl  +s2e  t/T2  (3) 

The  right-hand-side  of  the  above  equation  must  be  a  constant 
since  I  is  a  constant.  Thus  for  the  5  unknowns  £,Pi,P2.Vio,V2o,  there 
are  only  3  linearly  independent  equations,  namely,  (3)  and 
V10-R1I  =  s\,  V20-R2I  =  s2.  Therefore,  these  variables  cannot  be 
uniquely  determined. 

By  Claim  1,  we  know  that  the  circuit  parameters  can  not 
be  determined  with  unknown  initial  conditions  because  there  are 
not  sufficient  linearly  independent  equations  for  the  unknown 
variables. 


v(T0  +  Ti )  -  1 /(T0  +  2Ti )  =  yi  (Pi  -  l)pi  +  y2  (p2  -  1  )p2 

Since  pi,p2  are  given,  we  can  find  yi,  y2  from  the  above  two 
equations.  Since  si,s2,ai,a2  are  given,  we  can  compute  Ri,R2  from  (4). 
Then  Ci  =  iq/Pi,  C2  =  t2IR2  and  Vio  =  S1+P1L  v2o  =  s2+R2E  Finally, 
we  can  compute  parameter  E  from  (2). 

2.3.  Algorithms  for  parameter  identification 

For  completeness,  we  provide  the  algorithms  by  summarizing 
the  arguments  and  steps  from  the  previous  sections  and  using  the 
method  in  Ref.  [12], 

Let  the  load  current  be  /  for  0  <  t  <  To  and  f  for  t  <  To,  respec¬ 
tively.  The  terminal  voltage  v(t)  is  recorded  at  t  =  CT,  t  =  0+  and  for 
t  >  0.  Then  Ro  =  (v(0~)-v(0+))/L 

Algorithm  for  2nd-order  model.  Choose  T e  (0,To/4[  and  Ti  >  0. 

To  extract  the  7  parameters  E,Ri,R2,Ci,C2  and  vio,v2o,  we  will  need 
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the  measurement  of  the  terminal  voltage  v(t)  at  8  time  instants:  Let  the  roots  to  q 3  -  2u1q2  +  (u^  +  u2)q  +  (1/3  -  u-[U2)  —  0  be 

v(0+),v(T),v(2-r),  v(3r)tv(4r)Iv(7^)tv(To+T1)Iv(To+2T1).  qi,q2,<73.  Let 


Step  1:  Compute  s3  =  vw-Rih  s 2  =  v2o-R2I  and  t\  =  RiCi, 
T2  =  R2C2- 

For  fc  =  1,2, 3.4,  let  bu  =  v((k-l)T)-v(l<T). 

Compute 

L“2  J  L^J  “2J  L“4J 

Let  the  roots  to  the  second  order  equation  d2-U]d+u2  =  0  be 
di,d2.  Then 


Ul 

d2 

-b\ 

-1 

V 

u2 . 

b3 

—b2 

.v 

r  1  =  - 


ln(di)’ 


T2  =  “ 


ln(d2) 


Let 


*r 

1 

1 

-i 

V 

*2. 

di 

d2 

d2  . 

.  Then 


si 


x, 


di  -r 


S2  = 


*2 


d2  —  1 


Step  2:  Use  v(T^),  v(7o  +  Ti),  v(To  +  2Ti)  to  complete  the 
computation. 

Let  oq  =  e~r°/n ,  a2  =  e~r°/T2,  pi  =  e~r'/T1  ,p2  =  e-7'/’’2. 

Compute 


\d l] 

0 

1 

r 

-1 

<7i 

d2 

= 

1 

0 

1 

<72 

_  d3  . 

1 

1 

0 

<73 

Then 


ln(dj) 


J  =  1,2,3 


Let 


*1 

'  1 

1 

1  ' 

-1 

'V 

*2 

= 

di 

d2 

do 

^2 

*3. 

d2 

d3 

.^3. 

Then 


Step  2:  Let  aj  =  e  pj  =  e  7l/TJ,  j  =  1,2,3. 


yi 

pi  - 1  p2  - 1 

-1 

v(r+)-v(T0  +  h) 

Pi  -  1  P2  ~  1  P3  -  1 

Pl(Pl-l)  P2(P2-1). 

v(T0  +  h)-v(T0  +  2h) 

Form  Q  = 

Pl(Pl-l)  P2(P2-1)  P3(P3-1) 
.Pl(Pl-l)  P2(P2-1)  P3(P3-1)_ 

Then  Compute 


yi  -  s-iaq 

'"“Tut 


«2  = 


y2  -  S2a2 

yi 

=  Q-1 

v(r+)-v(T0  +  h) 

I-h  ’ 

y2 

v(T0  +  T,)-v(T0  +  2T1) 

Ly3  J 

Lv(T0  +  2T1)-v(T0  +  3r1)J 

Ci=r1/Ri,  C2  =  t2/R2 


Then 


Finally, 


■ho  =  si  +  l?i/,  v2q  =  s2  +  R2I, 


Rj  = 


yj-sJaj 

i-h  ’ 


Cj  =  Tj/Rj,  j  =  1,2,3 


and  £  =  v(0")  +  v10  +  v2o- 

For  3rd  or  higher  order  models,  the  algorithm  can  be  extended 
from  the  arguments  in  previous  sections  and  by  using  the  algebraic 
tool  developed  in  [12]. 

Algorithm  for  3rd-order  model.  Choose  Te(0,To/6]  and  Ti  >  0. 

To  extract  the  10  parameters  E,R\,R2<R3,C\,C2,Cs  and  V10.V20.V30,  we 
will  need  the  measurement  of  the  terminal  voltage  v(f)  at  11  time 
instants:  v(0+),v(/cT),k  =  1,2,. ..,6,  v(T+),  v(T0+kT^,  k  =  1,2,3. 

Step  1:  Compute  Sj  =  Vp—Rjl  and  tj  =  RjCj,  j  =  1,2,3  from 
v(0+),v(kT),fc  =  1,...,6. 

For  k  =  1,...,6,  let  £>/<  =  v((fc-l)T)-v(fcT). 

Compute 


Ul 

b3 

-b2 

K 

-1 

b4 

«2 

= 

b4 

-b3 

b2 

b5 

U3_ 

,b5 

-b4 

b3 

_  b6  _ 

Finally, 

Vjo  =  Sj  +  Rjl,  j=  1,2,3, 

and 

E  =  v(0_)  +  v\q  +  v20  +  v30. 

2.4.  Validation  of  the  algorithms  and  sensitivity  issues 

The  algorithms  can  be  easily  implemented  with  Matlab.  Four 
Matlab  codes  are  included  in  the  Appendix  for  verification  of  the 
algorithms  for  a  second-order  model  and  a  third-order  model,  with 
specific  explanations  for  the  notation  and  steps.  For  each  model, 
one  Matlab  code  is  used  to  generate  and  plot  the  voltage  response 
at  specified  time  instants,  given  all  the  parameters  and  initial 
conditions.  Another  Matlab  code  uses  the  voltage  response  to 
reconstruct  all  the  parameters  and  initial  conditions  by  using  the 
algorithm. 
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Both  algorithms  were  validated  except  for  very  small  numerical 
errors.  Here  we  use  the  computational  results  for  the  3rd-order 
model  to  demonstrate  the  key  algebraic  steps.  There  are  three  3  by 
3  matrices  in  the  Algorithm  which  need  to  take  inverse.  Their  ei¬ 
genvalues  are  given  below: 


eig 


( 


b3 

b4 

b5 


{0.0486,  -0.0217,  -0.0009} 


eig 


1 

d, 

d? 


1 

di 

4 


{1.1960,0.6844,0.0964} 


eig(Q)  =  {-0.9864,-0.2682,-0.0392} 

The  original  parameters  and  the  reconstructed  are  compared  as 
follows 


«1 

Ri 

r3 

Ci 

c2 

c3 

!'io 

!'20 

^30 


Original 

0.060000000000000 

0.050000000000000 

0.040000000000000 

200.000000000000 

640.000000000000 

2800.000000000000 

0.010000000000000 

0.020000000000000 

-0.050000000000000 


Reconstructed 

0.060000000000630 

0.050000000000336 

0.039999999998466 

199.999999998905 

640.000000007105 

2800.000000131836 

0.010000000000088 

0.020000000000599 

-0.050000000002005 


The  average  relative  error  (defined  as  jl-reconstructed/origi- 
nal[)  is  1.98  x  10 

For  the  second-order  model,  the  average  relative  error  is  much 
smaller,  7.08  x  1CT13. 

For  a  real  voltage  response  recorded  in  the  experiment,  mea¬ 
surement  error  and  noises  are  unavoidable.  Hence  it  is  necessary  to 
examine  the  sensitivity  of  the  reconstructed  parameters  with 
respect  to  voltage  perturbation.  If  all  voltage  values  are  increased 
(or  decreased)  at  the  same  time,  the  reconstructed  parameters  do 
not  change  much.  The  parameters  are  more  sensitive  to  the  type  of 
perturbation  where  the  voltage  is  increased  at  one  instant  but 
decreased  the  next.  So  we  used  the  following  perturbation  pattern: 
v(0+)+<5,  v(T)-i5,v(2T)+<5,  ...,  v(T+)  ~8,v(T0 +  TX)  +  6, ... 

Fig.  3  shows  some  relationship  between  the  average  relative 
parameter  deviation  and  the  perturbation  <5.  Clearly,  the  third-order 
model  is  much  more  sensitive  than  the  second-order  model.  On  the 
curve  for  the  3rd-order  model,  there  is  a  non-smooth  point,  which 
occurs  at  5  =  1.56  x  1CT5.  At  this  perturbation  value,  the  recon¬ 
structed  Ci  is  a  complex  number,  which  is  useless.  For  the  2nd- 
order  model,  complex  capacitance  will  also  occur,  but  at  a  much 
larger  5:  5  =  5.03  x  1CT4. 

Due  to  these  sensitivity  issues,  when  the  algorithm  for  a  3rd- 
order  model  is  used  on  experimental  responses,  the  reconstructed 
circuit  parameters  usually  have  complex  numbers.  When  the  al¬ 
gorithm  for  the  2nd-order  model  is  applied,  the  circuit  parameters 
are  usually  real  numbers  and  we  can  always  adjust  the  time  interval 
T  and  Ti  to  produce  reasonable  parameters  with  small  root-mean- 
square  error  between  experimental  response  and  model  response. 

2.5.  Considerations  for  extracting  circuit  parameters  from 
experimental  responses 

If  the  voltage  response  v(t)  is  computed  from  an  ideal  circuit 
model,  then  all  the  circuit  parameters  can  be  reconstructed,  with 


small  numerical  errors,  from  several  points: 
v(0~),v(0+),v(kT),v(To+kTi),k  =  1,2,....  The  results  will  not  be 
affected  significantly  by  the  choices  of  T  and  Ti,  except  for  very 
small  numerical  errors.  However,  the  experimental  responses  have 
various  nonidealities,  such  as  measurement  errors  and  quantiza¬ 
tion  errors.  Furthermore,  the  current  is  not  a  strict  step  function 
and  has  jitters  or  small  ripple.  After  all,  the  circuit  model  is  only  an 
approximation  of  the  complex  chemical  processes  inside  a  battery. 

Because  of  all  these  nonidealities,  the  parameters  extracted 
from  the  experimental  responses  will  depend  on  T,  Ti  and  To.  For 
some  values  of  T  and  Ti,  the  algorithm  may  produce  negative  or 
even  complex  numbers  for  the  resistance  and  capacitance.  This 
kind  of  situation  would  be  encountered  more  frequently  when 
higher-order  models  are  extracted.  Thus  we  need  to  vary  T  and  Ti 
within  an  admissible  range  and  choose  the  resulting  circuit  pa¬ 
rameters  which  produce  the  minimal  root-mean-square-error 
(RMSE)  between  the  experimental  response  and  the  response 
reconstructed  from  the  model. 

For  simplicity,  we  may  set  Ti  =  T and  use  one  dimensional  sweep 
to  find  the  best  value  T. 

Another  consideration  is  the  choice  of  the  difference  between 
load  currents  I  and  h.  Since  the  parameters  depend  on  the  dis¬ 
charging  current,  the  difference  should  be  kept  as  small  as  possible. 
However,  with  very  small  differences,  the  computed  parameters 
may  be  too  sensitive  to  noises  and  measurement  errors  and  the 
results  would  vary  widely  from  test  to  test.  In  our  experiment, 
different  values  of  l-l i  were  tested  and  a  suitable  one  was  chosen. 

The  total  test  time  and  the  switch  time  To  depend  on  what  kind 
of  dynamic  properties  are  evaluated  and  the  capacity  of  the  battery. 
For  brief  transient  properties,  a  few  seconds  or  minutes  of  response 
may  be  sufficient;  for  long  time  discharging/charging  under  steady 
load  current,  hours  of  test  time  may  be  required. 

3.  Experiment  and  computational  results  on  a  lead-acid 
battery 

The  experiment  was  conducted  to  demonstrate  the  use  of  our 
algorithm  in  identifying  circuit  parameters  under  non-zero  initial 
condition  of  the  capacitor  voltage.  It  was  not  intended  for  the 
complete  evaluation  of  a  battery  under  all  operating  conditions. 
Thus  the  tests  were  only  run  unde  some  particular  operating  con¬ 
ditions,  such  as  a  specific  discharging  current,  state  of  charge  and 
temperature.  Although  the  parameters  under  different  operating 
conditions  can  be  used  together  to  derive  more  complex  nonlinear 
models,  this  is  not  in  the  scope  of  this  paper. 
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The  battery  used  for  the  experiment  was  a  lead  acid  battery 
rated  6V  and  12  Ah.  The  2nd-order  Thevenin’s  equivalent  circuit 
model  was  considered.  We  examined  how  the  circuit  parameters 
depended  on  the  discharging  current  under  a  certain  state  of 
charge. 

A  group  of  tests  were  conducted  under  room  temperature 
(25  °C).  For  each  test,  the  current  load  was  programmed  to  draw  a 
step  current  which  was  initially  /  and  then  switched  to  £  =  /  +  0.2A 
at  To.  The  value  of  /  varies  from  0.6  A  to  1.8  A.  Each  test  lasted  about 
2  min.  The  rest  time  between  two  subsequent  tests  was  between  1 
and  5  min.  In  earlier  similar  experiment  for  the  work  [12],  the  rest 
time  was  about  30  min. 

The  terminal  voltage  responses  under  the  step  current  were 
recorded  via  a  16-digit  data  acquisition  (DAQ)  device.  The  input 
voltage  range  was  -10  V  to  10  V.  Thus  the  resolution  was  20/ 
216  =  3.0518~4  V.  No  digital  or  analog  filter  was  used  to  process  the 
data/signals. 

With  the  test  data,  we  computed  parameters  for  a  2nd-order 
model  under  different  discharging  current. 

Fig.  4  shows  the  terminal  voltage  v(t)  for  a  step  discharging 
current  with  I  =  0.6A  for  te( 0,75)  and  h  =  0.8A  for  t  >  75  s.  The 
experimental  response  is  plotted  in  pink  (light  colored)  curve.  The 
black  curve  is  the  response  reconstructed  from  the  circuit  model 
whose  parameters  were  extracted  from  the  experimental  response 
by  our  proposed  algorithm. 

For  simplicity,  we  chose  the  same  value  for  T  and  Tp  The  best  T 
was  obtained  via  a  one  dimensional  sweep  within  a  certain  range 
so  that  the  RMSE  is  minimal. 

The  computational  result  gave  T  =  Ti  =  3.69  s,  RMSE  =  0.0008  V. 
The  extracted  parameters  and  initial  conditions  are  given  below: 

Rq  —  0.02670,  Rj  -  0.02580,  R2  —  0.04120 
C\  =  53.79F,  C2  =  474. 3F 

i/10  =  -0.001  IV,  v20  =  -0.01  57V,  E  =  6.15481/ 

If  zero  initial  conditions  were  assumed,  we  obtained  different 
parameters  (except  for  R0),  as  given  below, 

R,  =  0.02890,  R2  =  0.06680 


C1  =  50.9F,  C2  =  308. IF 

Fig.  5  shows  the  experimental  response  for  I  =  1.8A  and  £  =  2 A 
with  switching  time  To  =  59.09  s.  The  computational  result  gives 
T  =  T]  =  4.905  s,  RMSE  =  0.0017.  The  larger  RMSE  is  due  to  larger 
discharging  current  and  larger  voltage  drop.  The  parameters  and 
initial  conditions  are  given  below: 

R0  =  0.02750,  R,  =  0.01660,  R2  =  0.01330, 

C,  =  94.9734F,  C2  =  854.71 53F 

i/10  =  -0.0149V,  v2q  =  -0.0227V,  £  =  6.1107V 

When  zero  initial  conditions  are  assumed,  we  obtained  the 
following  parameters, 

Rj  =  0.02460,  R2  =  0.02620 
Ci=61.44F,  C2  =  428. 9F 


Fig.  4.  Discharging  responses:  /  =  0.6/1,  /,  =  0.8/1. 


Fig.  6.  Parameters  vs  discharging  current. 
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Fig.  7.  Parameters  vs  discharging  current,  0  initial  condition  assumed. 


We  obtained  experimental  responses  for  I  =  0.6, 0.8, 1.0, 1.2, 
1.4, 1.8/4  and  h  =  I+0.2A.  Circuit  parameters  were  computed  for  each 
current.  Fig.  6  plots  the  computed  circuit  parameters  vs  the  dis¬ 
charging  current  I.  There  are  visible  changes  in  i?2-  Other  parame¬ 


ters  remain  nearly  flat.  The  curves  show  a  clear  pattern  of  the 
parameters  and  good  consistency  among  the  group  of  tests.  The 
curves  are  similar  to  those  in  Fig.  11  in  [12],  which  were  obtained  for 
a  different  lead-acid  battery. 

If  zero  initial  conditions  were  assumed,  different  parameters 
would  be  produced.  Fig.  7  show  the  result  when  zero  initial  con¬ 
ditions  are  assumed  for  the  computation. 

There  are  more  visible  changes  on  the  curves  and  the  patterns  of 
the  curves  are  not  as  clear  as  those  in  Fig.  6.  These  numbers  may 
lead  to  more  complex  relationship  and  functions  for  the  depen¬ 
dence  of  the  parameters  on  the  current. 


4.  Conclusions 

This  paper  derived  an  algebraic  method  for  identifying  param¬ 
eters  for  circuit  models  of  batteries  with  unknown  non-zero  initial 
conditions.  It  was  shown  that  such  a  parameter  identification 
problem  can  not  be  solved  with  a  constant  load  current  and  that  a 
step  load  current  is  necessary.  Since  the  new  method  does  not 
require  zero  initial  condition,  the  rest  time  between  two  subse¬ 
quent  tests  can  be  substantially  reduced,  thus  accelerating  the 
testing  cycles  for  battery  evaluation. 


A.  Matlab  codes 


Al.  Matlab  codes  for  a  2nd-order  model 

A.1.1.  Matlab  code  for  generating  voltage  response  under  step  current 

"/.Circuit  parameters.  R1C1,R2C2  in  ascending  order 
E=6 ; 

R0=0 . 02 ; 

Cl=200 ; R1=0 . 06 ; 


C2=2800 ; R2=0 . 04 ; 

“/.Initial  conditions 
vl0=0.02;v20=-0.05; 

“/.Values  for  step  current  and  time  intervals 
1=2;  11=2.2;  T0=200;  “/.I  for  t<T0,  II  for  t>T0 
T=30;T1=30; 


“/.exponential  constants:  d_i , alpha_i ,p_i ,  i=l,2 
dml=exp(-T/Rl/Cl) ;af l=exp(-T0/Rl/Cl) ;pl=exp(-Tl/Rl/Cl) ; 
dm2=exp (-T/R2/C2) ; af 2=exp(-T0/R2/C2) ;p2=exp(-Tl/R2/C2) ; 
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v0=E-vl0-v20-R0*I ;  "/.stands  for  v(0+) 
s=E- (R0+R1+R2) *1 ; 


vl=s- (vlO-Rl*I) *dml- (v20-R2*I) *dm2 ;  7.  for  v(T) 
v2=s- (vlO-Rl*I) *dml~ 2- (v20-R2*I) *dm2~2 
v3=s- (vlO-Rl*I) *dml‘3- (v20-R2*I) *dm2~3 
v4=s-  (vlO-Rl*I) *dml~4-  (v20-R2*I)  *dm2~4 ;  '/.for  v(4T) 


vT0=s- (vlO-Rl*I) *af 1- (v20-R2*I) *af 2 ;  '/.for  v(TO-) 

zl=(vlO-Rl*I) *af 1+R1* (I-Il) ; 

z2=(v20-R2*I) *af 2+R2* (I-Il) ; 

sl=E- (R0+R1+R2) *11 ; 

vT00=sl-zl-z2;  '/.for  v(T0+) 

vT01=sl-zl*pl-z2*p2 ;  "/.for  v(T0+Tl) 

vT02=sl-zl*pl~2-z2*p2~2 ;  '/.for  v(T0+2Tl) 

tt=  [0  0  T  T*2  T*3  T*4  TO  TO  T0+T1  T0+T1*2] ; 

vv= [E-vl0-v20  vO  vl  v2  v3  v4  vTO  vTOO  vTOl  vT02] ; 

'/, plot  (tt , vv)  '/of or  plotting  the  voltage  response 

'/.adding  perturbations 
z=0  '/.for  no  perturbation 
vO=vO+z ; vl=vl-z ; 


v2=v2+z ; v3=v3-z ; 
v4=v4+z ; vT00=vT00-z ; 
vT01=vT01+z ; vT02=vT02-z ; 
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A.1.2.  Matlab  code  for  reconstructing  parameters  and  initial  condition 

*/.  Given  values  of  voltage  response  at  specified  time  instants, 
l  step  discharging  current  I,  II,  with  stepping  time  TO, 

°/0  T  and  T1  are  time  intervals  before/after  TO 

°/0Step  1:  Computing  dl,d2,sl,s2 

bl=vO-vl;  b2=vl-v2;  °/0vO  for  v(0+) 

b3=v2-v3;  b4=v3-v4;  %vn  for  v(nT) 

u=inv([b2  -bl;b3  -b2] ) * [b3 ;b4] ; 

d=roots([l  -u(l)  u (2) ] ) ; 
dl=min(d) ;d2=max(d) ; 

x=inv([l  1;  dl  d2] ) *  [bl ; b2] ; 
sl=x(l) / (dl-1) ;s2=x(2)/(d2-l) ; 

°/0Step  2:  Reconstructing  parameters 
pl=dl~ (Tl/T) ;  p2=d2~ (Tl/T) ; 
af l=dl~ (TO/T) ;  af 2=d2~ (TO/T) ; 

A=  [pl-1  p2-l ;pl* (pl-1)  p2* (p2— 1 ) ] ; 

B= [vTOO-vTOl ; vT01-vT02] ; 
y=inv(A)*B; 

1  Reconstruct  R1 ,  R2 

Rel=(y (1) -sl*af 1) / (I-Il) ; 

Re2=(y(2)-s2*af2)/(I-Il) ; 

“/Reconstruct  Cl,  C2 
Cel=-T/Rel/log(dl) ; 

Ce2=-T/Re2/log(d2)  ; 
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“/.Reconstruct  vlO,  v20,  E; 
velO=Rel*I-x(l)/(l-dl) ; 
ve20=Re2*I-x(2)/(l-d2) ; 

E=v0+R0*I+vel0+ve20 ; 

“/.Relative  error  between  given  parameters  and  reconstructed  one 
“/»  RO  excluded  since  it  only  depends  on  v(O-)  and  v(0+) 

ee= [1-velO/vlO ; l-ve20/v20 ; 1-Rel/Rl ; 1-Re2/R2 ; 1-Cel/Cl ; 1-Ce2/C2 ; l-E/6]  ; 
“/.average  of  relative  error 
avgee=sum(abs (ee) ) /7 ; 

A.2.  Matlab  codes  for  a  3rd-order  model 

A2A.  Generating  voltage  response  under  step  current 

“/.Circuit  parameters.  R1C1  ,R2C2,R3C3  in  ascending  order 
E=6 ; 

R0=0 . 02 ; 

C1=200;R1=0.06; 

C2=640 ; R2=0 . 05 ; 

C3=2800 ; R3=0 . 04 ; 

“/.Initial  conditions 

vl0=0 . 01 ; v20=0 . 02 ; v30=-0 . 05 ; 

“/.Values  for  step  current  and  time  intervals 

1=2;  11=2. 2;  T0=200;  °/„I  for  t<T0,  II  for  t>T0 
T=30 ; Tl=30 ; 


“/.exponential  constants:  d_i,alpha_i,p_i,  i=l,2,3 

dml=exp(-T/Rl/Cl) ;af l=exp(-T0/Rl/Cl) ;pl=exp(-Tl/Rl/Cl) 
dm2=exp(-T/R2/C2) ; af 2=exp (-T0/R2/C2) ;p2=exp(-Tl/R2/C2) 
dm3=exp(-T/R3/C3) ; af 3=exp(-T0/R3/C3) ; p3=exp (-T1/R3/C3) 


“/.Computing  voltage  response  at  0+,T,2T,  .  .  .  ,T0,T0+T1  .T0+2T1 ,  .  .  . 
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v=zeros (1 , 12) ; 
s=E- (R0+R1+R2+R3) *1 ; 

for  n=l:7  %v(n)  stands  for  v((n-l)T) 

v(n)=s- (vlO-Rl*I) *dml~ (n-l)-(v20-R2*I) *dm2~ (n-1)- (v30-R3*I) *dm3~ (n-1) ; 
end 

v(8)=s-(vlO-Rl*I)*af l-(v20-R2*I)*af2-(v30-R3*I)*af3;  %v(8)  for  v(TO-); 

zl=(vlO-Rl*I)*af 1+R1*(I-I1) ; 

z2=(v20-R2*I) *af 2+R2* (I-Il) ; 

z3=(v30-R3*I) *af 3+R3* (I-Il) ; 

sl=E- (R0+R1+R2+R3) *11; 

v(9)=sl-zl-z2-z3;  %v(T0+) 

v(10)=sl-zl*pl-z2*p2-z3*p3;  %v(T0+Tl) 

v(ll)=sl-zl*pl~2-z2*p2~2-z3*p3~2;  %v(T0+2Tl) 

v(12)=sl-zl*pl~3-z2*p2~3-z3*p3~3;  %v(T0+3Tl) 

7.  Plotting  the  voltage  response 

tt=  [0  0  T  T*2  T*3  T*4  T*5  T*6  TO  TO  T0+T1  T0+T1*2  T0+T1*3] ; 
vv=  [E-vl0-v20-v30  v]  ;  °/0v(0-)=E-vl0-v20-v30 

%plot (tt , vv, ’ * ’ )  %  The  voltage  response  can  be  plotted 

%Add  some  perturbation  to  the  voltage  response 
ve=  [1  -11-11-11-11-11  -1]; 
delta=0;  ZChoose  a  small  number  for  delta 
v=v+ve*delta; 

%  The  vector  v  contains  all  the  voltage  values  at  11  specified  points 
7.  They  will  be  used  to  reproduce  the  circuit  parameters  along  with 
7„  I,  II,  TO,  T,  T1 

A.2.2.  Matlab  code  for  reconstructing  parameters  and  initial  condition 

7.  Given  values  of  voltage  response  at  specified  time  instants, 

7.  step  dischaging  current  I,  II,  with  stepping  time  TO, 

7o  T  and  T1  are  time  intervals  before/after  TO 
7«  Three  pairs  of  parallel  RC 
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“/.Step  1:  Computing  dl  ,d2  ,d3 ,  si ,  s2 ,  s3 

bl=v(l)-v(2) ;  b2=v(2)-v(3) ;  70v(l)  stands  for  v(0+) 

b3=v(3)-v(4) ;  b4=v(4)-v(5) ;  °/0v(n)  for  v((n-l)T); 

b5=v(5)-v(6) ;  b6=v(6)-v(7) ; 

u=inv([b3  -b2  bl;b4  -b3  b2;b5  -b4  b3] )*  [b4;b5;b6] ; 

q=roots([l  -2*u(l)  u(l)~2+u(2)  u(3)-u(l)*u(2)] ) ; 
d=inv(  [0  1  1 ; 1  0  1 ; 1  1  0])*q; 

d=sort(d);  “/.line  up  dl  d2  d3  in  ascending  order 
dl=d(l) ;d2=d(2) ;d3=d(3) ; 

x=inv(  [1  1  1;  dl  d2  d3;dl~2  d2~2  d3~2] ) *  [bl ;b2 ;b3] ; 
sl=x(l)/(dl-l) ;s2=x(2)/(d2-l) ;s3=x(3)/(d3-l) ; 

“/.Step  2:  Computing  circuit  parameter  and  initial  conditions 
pl=dl~ (Tl/T) ;  p2=d2~ (Tl/T) ;  p3=d3~ (Tl/T) ; 
af l=dl~ (TO/T) ;  af 2=d2~ (TO/T) ;  af 3=d3~ (TO/T) ; 

A=  [pl-l  p2-l  p3-l; 

pl*(pl-l)  p2*(p2-l)  p3*(p3-l); 
pl~2*(pl-l)  p2~2*(p2-l)  p3"2* (p3— 1 ) ] ; 

B=[v(9)-v(10) ;v(10)-v(ll) ;v(ll)-v(12)] ; 

°/»v(9)  for  v(T0+)  ;  v(10)  for  v(T0+Tl) 
y=inv(A)*B; 

°/«Re  construct  R1 ,  R2 ,  R3 

Rel=(y(l)-sl*af 1) / (I-Il) ; 

Re2=(y(2)-s2*af2)/(I-Il) ; 

Re3=(y(3)-s3*af3)/(I-Il) ; 

’/.  Reconstruct  Cl,  C2,  C3 

Cel=-T/Rel/log(dl) ;  ”/.T/log(dl)  for  taul 
Ce2=-T/Re2/log(d2)  ; 

Ce3=-T/Re3/log(d3) ; 


940 


L  Devarakonda,  T.  Hu  /  Journal  of  Power  Sources  268  (2014)  928-940 


’/.Reconstruct  initial  conditions  vlO,  v20,  v30  and  E 
velO=Rel*I-x(l)/ (1-dl) ; 
ve20=Re2*I-x(2)/(l-d2) ; 
ve30=Re3*I-x(3)/ (l-d3)  ; 
E=v(l)+R0*I+vel0+ve20+ve30; 


/.Relative  error  between  given  parameters  and  reconstructed  one 
°/„  R0  excluded  since  it  only  depends  on  v(O-)  and  v(0+) 


ee=  [1-velO/vlO ; l-ve20/v20 ; 1- 
1-Rel/Rl ; 1-Re2/R2 ; 1-Re3/R3; 
1-Cel/Cl ; 1-Ce2/C2 ; 1-Ce3/C3] 
’/.average  of  relative  deviation 
avgee=sum(abs (ee) ) / 10 ; 
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